package fem2.material;

import fem2.Element;
import fem2.element.StructuralElement;
import fem2.enu.State;
import inf.jlinalg.Array2DMatrix;
import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;

public class StVenantKirchhoffMM extends ElasticMaterialModel {

	private IMatrix De;

	public StVenantKirchhoffMM(double E, double nu, double thickness, double rho, State ss) {
		super(E, nu, thickness, rho, ss);
		/*
		 * set the constitutive matrix
		 */
		switch (state) {
		case PLANE_STRAIN:
			De = createPlaneStrainElasticityMatrix(E, nu);
			break;
		case PLANE_STRESS:
			De = createPlaneStressElasticityMatrix(E, nu);
			break;
		case THREE_DIMENSIONAL:
			De = create3dElasticityMatrix(E, nu);
			break;
		case ONE_DIMENSIONAL:
			De = new Array2DMatrix(1, 1);
			De.set(0, 0, E);
			break;
		default:
			throw new Error("invalid stress state");
		}
	}

	@Override
	public IMatrix getConstitutiveMatrix(Element e, int k) {
		return De;
	}

	@Override
	public void Initialize(Element e) {
	}

	@Override
	public void InitializeSolutionStep(Element e) {
		/**
		 * Calculate first stress for the element
		 */
		StructuralElement se = (StructuralElement) e;
		for (int i = 0; i < se.getNumIntegrationPoints(); i++) {
			IVector sigma = se.getStress(i);
			double[] xi = se.getIntegrationPoint(i);
			computeStress(se, xi, sigma);
		}
	}

	@Override
	public void InitalizeNonlinearIteration(Element e) {
	}

	@Override
	public void FinalizeNonlinearIteration(Element e) {
		// System.out.println("material finalize called");
		/*
		 * Calculate material response
		 */
		StructuralElement se = (StructuralElement) e;
		for (int i = 0; i < se.getNumIntegrationPoints(); i++) {
			IVector sigma = se.getStress(i);
			double[] xi = se.getIntegrationPoint(i);
			computeStress(se, xi, sigma);
		}
	}

	@Override
	public void FinalizeSolutionStep(Element e) {
	}

	@Override
	public void CalculateOutputData(Element e) {
	}

	/**
	 * Computes the stress at Gauss point k by standard operation sigma =
	 * C*epsilon where epsilon is total strain
	 * 
	 * @param se
	 *            element to provide strain routine
	 * @param sigma
	 *            stress need to be computed
	 * @param xi
	 *            Gauss point
	 */
	private void computeStress(StructuralElement se, double[] xi, IVector sigma) {
		IVector e = se.computeStrain(xi);
		BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, De, e, 0.0, sigma);
	}

	@Override
	public double computeEnergyAt(IVector e) {
		IVector s = (IVector) e.clone();
		BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, De, e, 0.0, s);
		return computeEnergyAt(e, s);
	}

	@Override
	public double computeEnergyAt(IVector e, IVector s) {
		return BLAM.dot(e, s) / 2;
	}
}
